/*
 * Matrix.h
 *
 * Created 5/10/2009 By Johnny Huynh
 *
 * Version 00.00.03 6/28/2009
 *
 * Copyright Information:
 * All content copyright © 2009 Johnny Huynh. All rights reserved.
 */
 
 // Matrix.h contains all the matrix functions

 #ifndef MATRIX_H
 #define MATRIX_H
 
 #include "macro.h"
 #include "Vector3.h"
 #include "OpenGL_Headers.h"
 #include <math.h>


 ///////////////////////////////// MATRIX 2 x 2 //////////////////////////////////////
 
 //#ifndef Matrix2f
 //#define Matrix2f Matrix2<GLfloat>
 //#endif // Matrix2f
 
 //#ifndef Matrix2d
 //#define Matrix2d Matrix2<double>
 //#endif // Matrix2d
 
 template <typename TYPENAME> class Matrix2;
 
 typedef Matrix2<GLfloat> Matrix2f;
 typedef Matrix2<GLdouble> Matrix2d;
 
 // The matrix based on the row and column
 //
 // | Xx Yx |
 // | Xy Yy |
 
 template <typename TYPENAME>
 class Matrix2
 {
 // Data Members
 public:
    TYPENAME Xx, Xy, Yx, Yy;    // contained in order of OpenGL notation (Xx, Xy, Yx, Yy)
 
 // Local Functions
    Matrix2( const TYPENAME& mXx = ZERO, const TYPENAME& mYx = ZERO,
             const TYPENAME& mXy = ZERO, const TYPENAME& mYy = ZERO );
    ~Matrix2();
    
 // Friend Functions
    template <typename TYPENAME> friend inline TYPENAME det( const Matrix2<TYPENAME>& m );
    
 };



 /** LOCAL FUNCTIONS **/
 
 /**
  * Constructor
  */
 template <typename TYPENAME>
 Matrix2<TYPENAME>::Matrix2( const TYPENAME& mXx = ZERO, const TYPENAME& mYx = ZERO,
                             const TYPENAME& mXy = ZERO, const TYPENAME& mYy = ZERO )
                            : Xx( mXx ), Xy( mXy ), Yx( mYx ), Yy( mYy )
 {
 
 }
 
 /**
  * Destructor
  */
 template <typename TYPENAME>
 Matrix2<TYPENAME>::~Matrix2()
 {
 
 }

 /**
  * det() returns the determinant of this matrix.
  *
  * @param (Matrix2<TYPENAME>&)m
  * @return TYPENAME
  */
 template <typename TYPENAME>
 inline TYPENAME det( const Matrix2<TYPENAME>& m )
 {
    return (m.Xx*m.Yy) - (m.Xy*m.Yx);
 }



 ///////////////////////////////// MATRIX 3 x 3 //////////////////////////////////////
 
 //#ifndef Matrix3f
 //#define Matrix3f Matrix3<GLfloat>
 //#endif // Matrix3f
 
 //#ifndef Matrix3d
 //#define Matrix3d Matrix3<double>
 //#endif // Matrix3d
 
 template <typename TYPENAME> class Matrix3;
 
 typedef Matrix3<GLfloat> Matrix3f;
 typedef Matrix3<double> Matrix3d;
 
 // The matrix based on the row and column
 //
 // | Xx Yx Zx |
 // | Xy Yy Zy |
 // | Xz Yz Zz |
 
 template <typename TYPENAME>
 class Matrix3
 {
 // Data Members
 public:
    TYPENAME Xx, Xy, Xz,
             Yx, Yy, Yz,
             Zx, Zy, Zz;    // stored in order of OpenGL notation (Xx, Xy, Xz, Yx, Yy, Yz, Zx, Zy, Zz)
 
 // Local Functions
    Matrix3( const TYPENAME& mXx = ZERO, const TYPENAME& mYx = ZERO, const TYPENAME& mZx = ZERO,
             const TYPENAME& mXy = ZERO, const TYPENAME& mYy = ZERO, const TYPENAME& mZy = ZERO,
             const TYPENAME& mXz = ZERO, const TYPENAME& mYz = ZERO, const TYPENAME& mZz = ZERO );
    Matrix3( const Matrix3<TYPENAME>& m );
    
 #ifdef VECTOR3_H
    Matrix3( const Vector3<TYPENAME>& u, const Vector3<TYPENAME>& v, const Vector3<TYPENAME>& t );
 #endif // VECTOR3_H
    
    ~Matrix3();
    inline Matrix3<TYPENAME>& operator=( const Matrix3<TYPENAME>& m );
    inline Matrix3<TYPENAME>& operator+=( const Matrix3<TYPENAME>& m );
    inline Matrix3<TYPENAME>& operator-=( const Matrix3<TYPENAME>& m );
    inline Matrix3<TYPENAME>& operator*=( const TYPENAME& s );
    inline Matrix3<TYPENAME>& operator*=( const Matrix3<TYPENAME>& m );
    inline Matrix3<TYPENAME>& operator/=( const TYPENAME& s );
    inline TYPENAME& operator[]( const size_t& i );
    inline const TYPENAME& operator[]( const size_t& i ) const;
    //inline TYPENAME& operator()( const GLint& i, const GLint& j );
    //inline const TYPENAME& operator()( const GLint& i, const GLint& j ) const;
    inline Matrix3<TYPENAME>& copy( const Matrix3<TYPENAME>& m );
    inline GLboolean isZeroMatrix() const;
    inline const TYPENAME * getArray() const;
    inline Matrix3<TYPENAME>& normalizeAxes();
    inline Matrix3<TYPENAME> getNormalizedAxes() const;
    inline Matrix3<TYPENAME>& transpose();
    inline Matrix3<TYPENAME> getTranspose() const;
    inline Matrix3<TYPENAME>& invert();
    inline Matrix3<TYPENAME> getInverse() const;
    
 // Friend Functions
    template <typename TYPENAME> friend inline bool operator==( const Matrix3<TYPENAME>& m, const Matrix3<TYPENAME>& n );
    template <typename TYPENAME> friend inline bool operator!=( const Matrix3<TYPENAME>& m, const Matrix3<TYPENAME>& n );
    template <typename TYPENAME> friend inline Matrix3<TYPENAME> operator+( const Matrix3<TYPENAME>& m, const Matrix3<TYPENAME>& n );
    template <typename TYPENAME> friend inline Matrix3<TYPENAME> operator-( const Matrix3<TYPENAME>& m, const Matrix3<TYPENAME>& n );
    template <typename TYPENAME> friend inline Matrix3<TYPENAME> operator-( const Matrix3<TYPENAME>& m );
    template <typename TYPENAME> friend inline Matrix3<TYPENAME> operator*( const TYPENAME& s, const Matrix3<TYPENAME>& m );
    template <typename TYPENAME> friend inline Matrix3<TYPENAME> operator*( const Matrix3<TYPENAME>& m, const TYPENAME& s );
    template <typename TYPENAME> friend inline Matrix3<TYPENAME> operator*( const Matrix3<TYPENAME>& m, const Matrix3<TYPENAME>& n );
    template <typename TYPENAME> friend inline Matrix3<TYPENAME> operator/( const Matrix3<TYPENAME>& m, const TYPENAME& s );
    template <typename TYPENAME> friend inline Matrix3<TYPENAME> inverse( const Matrix3<TYPENAME>& m );
    template <typename TYPENAME> friend inline Matrix3<TYPENAME> normalAxes( const Matrix3<TYPENAME>& m );
    template <typename TYPENAME> friend inline Matrix3<TYPENAME>& normalizeAxes( Matrix3<TYPENAME>& m );
    template <typename TYPENAME> friend inline Matrix3<TYPENAME> transpose( const Matrix3<TYPENAME>& m );
    template <typename TYPENAME> friend inline TYPENAME det( const Matrix3<TYPENAME>& m );
    template <typename TYPENAME> friend inline Matrix3<TYPENAME> adjoint( const Matrix3<TYPENAME>& m );
    template <typename TYPENAME> friend inline void print( const Matrix3<TYPENAME>& m );
    
    #ifdef VECTOR3_H
    template <typename TYPENAME> friend inline Vector3<TYPENAME> operator*( const Matrix3<TYPENAME>& m, const Vector3<TYPENAME>& u );
    template <typename TYPENAME> friend inline Vector3<TYPENAME> operator*( const Vector3<TYPENAME>& u, const Matrix3<TYPENAME>& m );
    #endif // VECTOR3_H
    
 };



 /** LOCAL FUNCTIONS **/
 
 /**
  * Constructor
  */
 template <typename TYPENAME>
 Matrix3<TYPENAME>::Matrix3( const TYPENAME& mXx= ZERO, const TYPENAME& mYx = ZERO, const TYPENAME& mZx = ZERO,
                             const TYPENAME& mXy= ZERO, const TYPENAME& mYy = ZERO, const TYPENAME& mZy = ZERO,
                             const TYPENAME& mXz= ZERO, const TYPENAME& mYz = ZERO, const TYPENAME& mZz = ZERO )
                           : Xx( mXx ), Xy( mXy ), Xz( mXz ),
                             Yx( mYx ), Yy( mYy ), Yz( mYz ),
                             Zx( mZx ), Zy( mZy ), Zz( mZz )
 {

 }
 
 /**
  * Alternative Constructor
  */
 template <typename TYPENAME>
 Matrix3<TYPENAME>::Matrix3( const Matrix3<TYPENAME>& m )
 {
    memcpy( this, &m, sizeof(Matrix3<TYPENAME>) );
 }
 
 #ifdef VECTOR3_H
 /**
  * Alternative Constructor
  *
  *     [u.x v.x t.x
  * M =  u.y v.y t.y
  *      u.z v.z t.z]
  */
 template <typename TYPENAME>
 Matrix3<TYPENAME>::Matrix3( const Vector3<TYPENAME>& u, const Vector3<TYPENAME>& v, const Vector3<TYPENAME>& t )
                           : Xx( u.x ), Xy( u.y ), Xz( u.z ),
                             Yx( v.x ), Yy( v.y ), Yz( v.z ),
                             Zx( t.x ), Zy( t.y ), Zz( t.z )
 {

 }
 #endif // VECTOR3_H
 
 /**
  * Destructor
  */
 template <typename TYPENAME>
 Matrix3<TYPENAME>::~Matrix3()
 {
 
 }
 
 /**
  * operator=() copies the values of the matrix m
  * and returns this matrix (referenced).
  * It returns a reference to itself, not matrix m.
  *
  * @param (const Matrix3<TYPENAME>&)m
  * @return Matrix3<TYPENAME>&
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME>& Matrix3<TYPENAME>::operator=( const Matrix3<TYPENAME>& m )
 {  
    memcpy( this, &m, sizeof(Matrix3<TYPENAME>) );
    return *this;
 }

 /**
  * operator+=() add the values of the matrix m
  * to this matrix and returns this matrix (referenced).
  *
  * @param (Matrix3<TYPENAME>)m
  * @return Matrix3<TYPENAME>&
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME>& Matrix3<TYPENAME>::operator+=( const Matrix3<TYPENAME>& m )
 {
    Xx += m.Xx;
    Xy += m.Xy;
    Xz += m.Xz;
    Yx += m.Yx;
    Yy += m.Yy;
    Yz += m.Yz;
    Zx += m.Zx;
    Zy += m.Zy;
    Zz += m.Zz;
    return *this;
 }
 
 /**
  * operator-=() subtracts the values of the matrix m
  * from this matrix and returns this matrix (referenced).
  *
  * @param (Matrix3<TYPENAME>)m
  * @return Matrix3<TYPENAME>&
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME>& Matrix3<TYPENAME>::operator-=( const Matrix3<TYPENAME>& m )
 {
    Xx -= m.Xx;
    Xy -= m.Xy;
    Xz -= m.Xz;
    Yx -= m.Yx;
    Yy -= m.Yy;
    Yz -= m.Yz;
    Zx -= m.Zx;
    Zy -= m.Zy;
    Zz -= m.Zz;
    return *this;
 }
 
 /**
  * operator*=() returns the product of this matrix and the specified a scalar.
  * This matrix is also modified as a result of the multiplication.
  *
  * @param (const TYPENAME&) s
  * @return Matrix3<TYPENAME>&
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME>& Matrix3<TYPENAME>::operator*=( const TYPENAME& s )
 {
    Xx *= s;
    Yx *= s;
    Zx *= s;
    Xy *= s;
    Yy *= s;
    Zy *= s;
    Xz *= s;
    Yz *= s;
    Zz *= s;
    
    return *this;
 }
 
 /**
  * operator*=() returns the product between this matrix and the supplied matrix.
  * This matrix is also modified as a result of the multiplication.
  *
  * @param (const Matrix3<TYPENAME>&) m
  * @return Matrix3<TYPENAME>&
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME>& Matrix3<TYPENAME>::operator*=( const Matrix3<TYPENAME>& m )
 {
    TYPENAME a = (Xx * m.Xx) + (Yx * m.Xy) + (Zx * m.Xz); // Xx
    TYPENAME b = (Xx * m.Yx) + (Yx * m.Yy) + (Zx * m.Yz); // Yx
    TYPENAME c = (Xx * m.Zx) + (Yx * m.Zy) + (Zx * m.Zz); // Zx
    Xx = a;
    Yx = b;
    Zx = c;
    
    a = (Xy * m.Xx) + (Yy * m.Xy) + (Zy * m.Xz); // Xy
    b = (Xy * m.Yx) + (Yy * m.Yy) + (Zy * m.Yz); // Yy
    c = (Xy * m.Zx) + (Yy * m.Zy) + (Zy * m.Zz); // Zy
    Xy = a;
    Yy = b;
    Zy = c;
    
    a = (Xz * m.Xx) + (Yz * m.Xy) + (Zz * m.Xz); // Xz
    b = (Xz * m.Yx) + (Yz * m.Yy) + (Zz * m.Yz); // Yz
    c = (Xz * m.Zx) + (Yz * m.Zy) + (Zz * m.Zz); // Zz
    Xz = a;
    Yz = b;
    Zz = c;
    
    return *this;
 }
 
 /**
  * operator/=() returns the division of a matrix by a scalar. This
  * matrix is also modified by the division.
  *
  * @param (const Matrix3<TYPENAME>&) m, (const TYPENAME&) s
  * @return Matrix3<TYPENAME>&
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME>& Matrix3<TYPENAME>::operator/=( const TYPENAME& s )
 {
    Xx /= s;
    Yx /= s;
    Zx /= s;
    Xy /= s;
    Yy /= s;
    Zy /= s;
    Xz /= s;
    Yz /= s;
    Zz /= s;
    
    return *this;
 }
 
 /**
  * operator[]() returns the ith element of the matrix.
  * Note that 0 <= i <= 8.
  *
  * @param (const size_t&) i
  * @return TYPENAME&
  */
 template <typename TYPENAME>
 inline TYPENAME& Matrix3<TYPENAME>::operator[]( const size_t& i )
 {
    return reinterpret_cast<TYPENAME*>(this)[i];
 }
 
 /**
  * operator[]() returns the ith element of the matrix.
  * Note that 0 <= i <= 8.
  *
  * @param (const size_t&) i
  * @return const TYPENAME&
  */
 template <typename TYPENAME>
 inline const TYPENAME& Matrix3<TYPENAME>::operator[]( const size_t& i ) const
 {
    return reinterpret_cast<TYPENAME*>(this)[i];
 }

 /**
  * operator()() returns the element at the ith row and jth column of the matrix.
  * Note that 0 <= i, j <= 2.
  *
  * @param (const size_t&) i, (const size_t&) j
  * @return TYPENAME&
  */
 /*template <typename TYPENAME>
 inline TYPENAME& Matrix3<TYPENAME>::operator()( const size_t& i, const size_t& j )
 {
    return reinterpret_cast<TYPENAME*>(this)[i + (j*0x3)];
 }*/
 
 /**
  * operator()() returns the element at the ith row and jth column of the matrix.
  * Note that 0 <= i, j <= 2.
  *
  * @param (const size_t&) i, (const size_t&) j
  * @return const TYPENAME&
  */
 /*template <typename TYPENAME>
 inline const TYPENAME& Matrix3<TYPENAME>::operator()( const size_t& i, const size_t& j ) const
 {
    return reinterpret_cast<TYPENAME*>(this)[i + (j*0x3)];
 }*/

/**
  * copy() copies the values of the matrix m
  * and returns this matrix (referenced).
  * It returns a reference to itself, not matrix m.
  *
  * @param (Matrix3<TYPENAME>) u
  * @return Matrix3<TYPENAME>&
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME>& Matrix3<TYPENAME>::copy( const Matrix3<TYPENAME>& m )
 {
    memcpy( this, &m, sizeof(Matrix3<TYPENAME>) );
    return *this;
 }
 
 /**
  * isZeroMatrix() returns true if this matrix is the zero matrix;
  * otherwise, false is returned.
  *
  * @return GLboolean
  */
 template <typename TYPENAME>
 inline GLboolean Matrix3<TYPENAME>::isZeroMatrix() const
 {
    return Xx == ZERO && Xy == ZERO && Xz == ZERO
            && Yx == ZERO && Yy == ZERO && Yz == ZERO
            && Zx == ZERO && Zy == ZERO && Zz == ZERO;
 }

 /**
  * getArray() returns the OpenGL array representation of this matrix 
  * (i.e. [Xx, Xy, Xz, Yx, Yy, Yz, Zx, Zy, Zz]).
  *
  * @return const TYPENAME *
  */
 template <typename TYPENAME>
 inline const TYPENAME * Matrix3<TYPENAME>::getArray() const
 {
    return reinterpret_cast<const TYPENAME *>(this);
 }

 /**
  * normalizeAxes() converts this matrix into a matrix with unit vectors
  * along each column and returns a reference to this matrix.
  *
  * @return Matrix3<TYPENAME>&
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME>& Matrix3<TYPENAME>::normalizeAxes()
 {
    // normalize X-axis
    TYPENAME m( (Xx*Xx) + (Xy*Xy) + (Xz*Xz) );
    if ( m != ZERO )
    {
        m = sqrt( m );
        Xx /= m;
        Xy /= m;
        Xz /= m;
    }
    
    // normalize Y-axis
    m = (Yx*Yx) + (Yy*Yy) + (Yz*Yz);
    if ( m != ZERO )
    {
        m = sqrt( m );
        Yx /= m;
        Yy /= m;
        Yz /= m;
    }
    
    // normalize Z-axis
    m = (Zx*Zx) + (Zy*Zy) + (Zz*Zz);
    if ( m != ZERO )
    {
        m = sqrt( m );
        Zx /= m;
        Zy /= m;
        Zz /= m;
    }
    
    return *this;
 }
 
 /**
  * getNormalizedAxes() returns a copy of this matrix with its
  * axes normalized.
  * Notice this matrix is not altered.
  * Instead, a newly created matrix is modified.
  *
  * @return Matrix3<TYPENAME>
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME> Matrix3<TYPENAME>::getNormalizedAxes() const
 {
    Matrix3<TYPENAME> n( *this );
    
    // normalize X-axis
    TYPENAME m( (n.Xx*n.Xx) + (n.Xy*n.Xy) + (n.Xz*n.Xz) );
    if ( m != ZERO )
    {
        m = sqrt( m );
        n.Xx /= m;
        n.Xy /= m;
        n.Xz /= m;
    }
    
    // normalize Y-axis
    m = (n.Yx*n.Yx) + (n.Yy*n.Yy) + (n.Yz*n.Yz);
    if ( m != ZERO )
    {
        m = sqrt( m );
        n.Yx /= m;
        n.Yy /= m;
        n.Yz /= m;
    }
    
    // normalize Z-axis
    m = (n.Zx*n.Zx) + (n.Zy*n.Zy) + (n.Zz*n.Zz);
    if ( m != ZERO )
    {
        m = sqrt( m );
        n.Zx /= m;
        n.Zy /= m;
        n.Zz /= m;
    }
    
    return n;
 }
 
 /**
  * transpose() modifies this matrix to become the transpose
  * of this matrix.
  *
  * @return Matrix3<TYPENAME>&
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME>& Matrix3<TYPENAME>::transpose()
 {
    TYPENAME temp( Xy );
    Xy = Yx;
    Yx = temp;
    
    temp = Xz;
    Xz = Zx;
    Zx = temp;
    
    temp = Yz;
    Yz = Zy;
    Zy = temp;
    
    return *this;
 }
 
 /**
  * getTranspose() returns the tranpose version of this matrix.
  * Notice this matrix is not altered.
  * Instead, a new matrix is created.
  *
  * @return Matrix3<TYPENAME>
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME> Matrix3<TYPENAME>::getTranspose() const
 {
    return Matrix3<TYPENAME>( Xx, Xy, Xz,
                              Yx, Yy, Yz,
                              Zx, Zy, Zz );
 }

 /**
  * invert() modifies this matrix to become the inverse
  * of this matrix.
  * If the inverse does not exist, the matrix is set to
  * a zero matrix.
  *
  * @return Matrix3<TYPENAME>&
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME>& Matrix3<TYPENAME>::invert()
 {
    // Inverse = adjoint / determinant
    TYPENAME determinant( det( *this ) );
    if ( determinant )  // if determinant is not equal to 0
        *this = adjoint( *this ) / determinant;
    else
    {
        Xx = ZERO;
        Xy = ZERO;
        Xz = ZERO;
        Yx = ZERO;
        Yy = ZERO;
        Yz = ZERO;
        Zx = ZERO;
        Zy = ZERO;
        Zz = ZERO;
    }
    
    return *this;
 }
 
 /**
  * getInverse() returns an inverted copy of this matrix.
  * Notice this matrix is not altered.
  * Instead, a newly created matrix is modified.
  * If the inverse does not exist, the zero matrix is returned.
  *
  * @return Matrix3<TYPENAME>
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME> Matrix3<TYPENAME>::getInverse() const
 {
    TYPENAME determinant = det( *this );
    if ( determinant )  // if determinant is not equal to 0
        return adjoint( *this ) / determinant;
    else
        return Matrix3<TYPENAME>( ZERO, ZERO, ZERO,
                                  ZERO, ZERO, ZERO,
                                  ZERO, ZERO, ZERO );
 }
 
 /** FRIEND FUNCTIONS **/
 
 /**
  * operator==() returns true if the two specified matrices are equal;
  * otherwise, false is returned.
  *
  * @param (const Matrix3<TYPENAME>&) m
  * @param (const Matrix3<TYPENAME>&) n
  * @return bool
  */
 template <typename TYPENAME>
 inline bool operator==( const Matrix3<TYPENAME>& m, const Matrix3<TYPENAME>& n )
 {  
    return m.Xx == n.Xx && m.Yx == n.Yx && m.Zx == n.Zx
        && m.Xy == n.Xy && m.Yy == n.Yy && m.Zy == n.Zy
        && m.Xz == n.Xz && m.Yz == n.Yz && m.Zz == n.Zz;
 }
 
 /**
  * operator!=() returns true if the two specified matrices are not equal;
  * otherwise, false is returned.
  *
  * @param (const Matrix3<TYPENAME>&) m
  * @param (const Matrix3<TYPENAME>&) n
  * @return bool
  */
 template <typename TYPENAME>
 inline bool operator!=( const Matrix3<TYPENAME>& m, const Matrix3<TYPENAME>& n )
 {  
    return m.Xx != n.Xx || m.Yx != n.Yx || m.Zx != n.Zx
        || m.Xy != n.Xy || m.Yy != n.Yy || m.Zy != n.Zy
        || m.Xz != n.Xz || m.Yz != n.Yz || m.Zz != n.Zz;
 }
 
 /**
  * operator+() returns a vector constituting the addition 
  * between the two supplied matrices (i.e. m + n).
  *
  * @param (const Matrix3<TYPENAME>&) m, (const Matrix3<TYPENAME>&) n
  * @return Matrix3<TYPENAME>
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME> operator+( const Matrix3<TYPENAME>& m, const Matrix3<TYPENAME>& n )
 {
    return Matrix3<TYPENAME>( m.Xx + n.Xx, m.Yx + n.Yx, m.Zx + n.Zx,
                              m.Xy + n.Xy, m.Yy + n.Yy, m.Zy + n.Zy,
                              m.Xz + n.Xz, m.Yz + n.Yz, m.Zz + n.Zz );
 }

 /**
  * operator-() returns a matrix constituting the subtraction 
  * between the two supplied matrices (i.e. m - n).
  *
  * @param (const Matrix3<TYPENAME>&) m, (const Matrix3<TYPENAME>&) n
  * @return Matrix3<TYPENAME>
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME> operator-( const Matrix3<TYPENAME>& m, const Matrix3<TYPENAME>& n )
 {
    return Matrix3<TYPENAME>( m.Xx - n.Xx, m.Yx - n.Yx, m.Zx - n.Zx,
                              m.Xy - n.Xy, m.Yy - n.Yy, m.Zy - n.Zy,
                              m.Xz - n.Xz, m.Yz - n.Yz, m.Zz - n.Zz );
 }

 /**
  * operator-() returns a matrix with the values
  * of matrix m negated (i.e. -m).
  *
  * @param (const Matrix3<TYPENAME>&)m
  * @return Matrix3<TYPENAME>
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME> operator-( const Matrix3<TYPENAME>& m )
 {
    return Matrix3<TYPENAME>( -m.Xx, -m.Yx, -m.Zx,
                              -m.Xy, -m.Yy, -m.Zy,
                              -m.Xz, -m.Yz, -m.Zz );
 }

 /**
  * operator*() returns the product of a scalar and a matrix.
  *
  * @param (const TYPENAME&) s, (const Matrix3<TYPENAME>&) m
  * @return Matrix3<TYPENAME>
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME> operator*( const TYPENAME& s, const Matrix3<TYPENAME>& m )
 {
    return Matrix3<TYPENAME>( s*m.Xx, s*m.Yx, s*m.Zx,
                              s*m.Xy, s*m.Yy, s*m.Zy,
                              s*m.Xz, s*m.Yz, s*m.Zz );
 }

 /**
  * operator*() returns the product of a scalar and a matrix.
  *
  * @param (const Matrix3<TYPENAME>&) m, (const TYPENAME&) s
  * @return Matrix3<TYPENAME>
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME> operator*( const Matrix3<TYPENAME>& m, const TYPENAME& s )
 {
    return s*m;
 }

 /**
  * operator*() returns the product between the two supplied matrices.
  *
  * @param (const Matrix3<TYPENAME>&) m, (const Matrix3<TYPENAME>&) n
  * @return Matrix3<TYPENAME>
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME> operator*( const Matrix3<TYPENAME>& m, const Matrix3<TYPENAME>& n )
 {
    return Matrix3<TYPENAME>( (m.Xx * n.Xx) + (m.Yx * n.Xy) + (m.Zx * n.Xz),   // Xx
                              (m.Xx * n.Yx) + (m.Yx * n.Yy) + (m.Zx * n.Yz),   // Yx
                              (m.Xx * n.Zx) + (m.Yx * n.Zy) + (m.Zx * n.Zz),   // Zx
                              (m.Xy * n.Xx) + (m.Yy * n.Xy) + (m.Zy * n.Xz),   // Xy
                              (m.Xy * n.Yx) + (m.Yy * n.Yy) + (m.Zy * n.Yz),   // Yy
                              (m.Xy * n.Zx) + (m.Yy * n.Zy) + (m.Zy * n.Zz),   // Zy
                              (m.Xz * n.Xx) + (m.Yz * n.Xy) + (m.Zz * n.Xz),   // Xz
                              (m.Xz * n.Yx) + (m.Yz * n.Yy) + (m.Zz * n.Yz),   // Yz
                              (m.Xz * n.Zx) + (m.Yz * n.Zy) + (m.Zz * n.Zz) ); // Zz
 }

 /**
  * operator/() returns the division of a matrix by a scalar.
  *
  * @param (const Matrix3<TYPENAME>&) m, (const TYPENAME&) s
  * @return Matrix3<TYPENAME>
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME> operator/( const Matrix3<TYPENAME>& m, const TYPENAME& s )
 {
    return Matrix3<TYPENAME>( m.Xx/s, m.Yx/s, m.Zx/s,
                              m.Xy/s, m.Yy/s, m.Zy/s,
                              m.Xz/s, m.Yz/s, m.Zz/s );
 }

 /**
  * inverse() returns the inverse matrix of the specified matrix.
  * If the inverse does not exist, the zero matrix is returned.
  *
  * @return Matrix3<TYPENAME>
  */
 template <typename TYPENAME> 
 inline Matrix3<TYPENAME> inverse( const Matrix3<TYPENAME>& m )
 {
    TYPENAME determinant( det( m ) );
    if ( determinant )  // if determinant is not equal to 0
        return adjoint( m ) / determinant;
    else
        return Matrix3<TYPENAME>( ZERO, ZERO, ZERO,
                                  ZERO, ZERO, ZERO,
                                  ZERO, ZERO, ZERO );
 }
 
 /**
  * normalAxes() returns a form of the specified matrix with its
  * axes normalized.
  *
  * @param (const Matrix3<TYPENAME>&) m
  * @return Matrix3<TYPENAME>
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME> orthogonal( const Matrix3<TYPENAME>& m )
 {
    Matrix3<TYPENAME> n( m );
    
    // normalize X-axis
    TYPENAME mag( (n.Xx*n.Xx) + (n.Xy*n.Xy) + (n.Xz*n.Xz) );
    if ( mag != ZERO )
    {
        mag = sqrt( mag );
        n.Xx /= mag;
        n.Xy /= mag;
        n.Xz /= mag;
    }
    
    // normalize Y-axis
    mag = (n.Yx*n.Yx) + (n.Yy*n.Yy) + (n.Yz*n.Yz);
    if ( mag != ZERO )
    {
        mag = sqrt( mag );
        n.Yx /= mag;
        n.Yy /= mag;
        n.Yz /= mag;
    }
    
    // normalize Z-axis
    mag = (n.Zx*n.Zx) + (n.Zy*n.Zy) + (n.Zz*n.Zz);
    if ( mag != ZERO )
    {
        mag = sqrt( mag );
        n.Zx /= mag;
        n.Zy /= mag;
        n.Zz /= mag;
    }
    
    return n;
 }
 
 /**
  * normalizeAxes() converts this matrix into a matrix with unit vectors
  * along each column and returns a reference to this matrix.
  *
  * @param (Matrix3<TYPENAME>&) m
  * @return Matrix3<TYPENAME>&
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME>& normalizeAxes( Matrix3<TYPENAME>& m )
 {
    // normalize X-axis
    TYPENAME mag( (m.Xx*m.Xx) + (m.Xy*m.Xy) + (m.Xz*m.Xz) );
    if ( mag != ZERO )
    {
        mag = sqrt( mag );
        m.Xx /= mag;
        m.Xy /= mag;
        m.Xz /= mag;
    }
    
    // normalize Y-axis
    mag = (m.Yx*m.Yx) + (m.Yy*m.Yy) + (m.Yz*m.Yz);
    if ( mag != ZERO )
    {
        mag = sqrt( mag );
        m.Yx /= mag;
        m.Yy /= mag;
        m.Yz /= mag;
    }
    
    // normalize Z-axis
    mag = (m.Zx*m.Zx) + (m.Zy*m.Zy) + (m.Zz*m.Zz);
    if ( mag != ZERO )
    {
        mag = sqrt( mag );
        m.Zx /= mag;
        m.Zy /= mag;
        m.Zz /= mag;
    }
    
    return m;
 }
 
 /**
  * transpose() returns the tranpose of the specified matrix.
  * Notice the specified matrix is not altered.
  * Instead, a new matrix is created.
  *
  * @param (const Matrix3<TYPENAME>&) m
  * @return Matrix3<TYPENAME>
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME> transpose( const Matrix3<TYPENAME>& m )
 {
    return Matrix3<TYPENAME>( m.Xx, m.Xy, m.Xz,
                              m.Yx, m.Yy, m.Yz,
                              m.Zx, m.Zy, m.Zz );
 }
 
 /**
  * det() returns the determinant of this matrix.
  *
  * @param (Matrix3<TYPENAME>&) m
  * @return TYPENAME
  */
 template <typename TYPENAME>
 inline TYPENAME det( const Matrix3<TYPENAME>& m )
 {
    // m.Xx * det( m.Yy, m.Zy   - m.Yx * det( m.Xy, m.Zy    + m.Zx * det( m.Xy, m.Yy
    //             m.Yz, m.Zz )               m.Xz, m.Zz )                m.Xz, m.Yz )
    return (m.Xx*( (m.Yy*m.Zz) - (m.Zy*m.Yz) ))
         - (m.Yx*( (m.Xy*m.Zz) - (m.Zy*m.Xz) ))
         + (m.Zx*( (m.Xy*m.Yz) - (m.Yy*m.Xz) ));
 }
 
 /**
  * adjoint() returns the adjoint matrix of the specified matrix.
  *
  * @param (const Matrix3<TYPENAME>&) m
  * @return Matrix3<TYPENAME>
  */
 template <typename TYPENAME>
 inline Matrix3<TYPENAME> adjoint( const Matrix3<TYPENAME>& m )
 {
    /*TYPENAME detXx = det( Matrix2<TYPENAME>(m.Yy, m.Zy,
                                            m.Yz, m.Zz) );
    TYPENAME detYx = -det( Matrix2<TYPENAME>(m.Xy, m.Zy,
                                             m.Xz, m.Zz) );
    TYPENAME detZx = det( Matrix2<TYPENAME>(m.Xy, m.Yy,
                                            m.Xz, m.Yz) );
    TYPENAME detXy = -det( Matrix2<TYPENAME>(m.Yx, m.Zx,
                                             m.Yz, m.Zz) );
    TYPENAME detYy = det( Matrix2<TYPENAME>(m.Xx, m.Zx,
                                            m.Xz, m.Zz) );
    TYPENAME detZy = -det( Matrix2<TYPENAME>(m.Xx, m.Yx,
                                             m.Xz, m.Yz) );
    TYPENAME detXz = det( Matrix2<TYPENAME>(m.Yx, m.Zx,
                                            m.Yy, m.Zy) );
    TYPENAME detYz = -det( Matrix2<TYPENAME>(m.Xx, m.Zx,
                                             m.Xy, m.Zy) );
    TYPENAME detZz = det( Matrix2<TYPENAME>(m.Xx, m.Yx,
                                            m.Xy, m.Yy) );
    
    // transpose of matrix with the determinant at each entry
    return Matrix3<TYPENAME>( detXx, detXy, detXz,
                              detYx, detYy, detYz,
                              detZx, detZy, detZz );*/
    
    return Matrix3<TYPENAME>( (m.Yy*m.Zz) - (m.Zy*m.Yz), (m.Zx*m.Yz) - (m.Yx*m.Zz), (m.Yx*m.Zy) - (m.Zx*m.Yy),
                              (m.Zy*m.Xz) - (m.Xy*m.Zz), (m.Xx*m.Zz) - (m.Zx*m.Xz), (m.Zx*m.Xy) - (m.Xx*m.Zy),
                              (m.Xy*m.Yz) - (m.Yy*m.Xz), (m.Yx*m.Xz) - (m.Xx*m.Yz), (m.Xx*m.Yy) - (m.Yx*m.Xy) );
 }
 
 /**
  * print() prints out information about the specified matrix.
  *
  * @param (const Matrix3<TYPENAME>&) m
  */
 template <typename TYPENAME>
 inline void print( const Matrix3<TYPENAME>& m )
 {
    printf("[%1f %2f %3f\n %4f %5f %6f\n %7f %8f %9f]\n", m.Xx, m.Yx, m.Zx, 
                                                      m.Xy, m.Yy, m.Zy, 
                                                      m.Xz, m.Yz, m.Zz );
 }

 #ifdef VECTOR3_H
 
 /**
  * operator*() returns a column vector by multiplying the specified matrix by a column vector.
  *
  * [Xx Yx Zx      [x
  *  Xy Yy Zy   *   y
  *  Xz Yz Zz]      z]
  *
  * @param (const Matrix3<TYPENAME>&) m - matrix
  * @param (const Vector3<TYPENAME>&) u - column vector
  * @return Vector3<TYPENAME>
  */
 template <typename TYPENAME>
 inline Vector3<TYPENAME> operator*( const Matrix3<TYPENAME>& m, const Vector3<TYPENAME>& u )
 {
    return Vector3<TYPENAME>( (u.x*m.Xx) + (u.y*m.Yx) + (u.z*m.Zx), 
                              (u.x*m.Xy) + (u.y*m.Yy) + (u.z*m.Zy), 
                              (u.x*m.Xz) + (u.y*m.Yz) + (u.z*m.Zz) );
 }
 
 /**
  * operator*() returns a row vector by multiplying the specified row vector by a matrix.
  *
  *            [Xx Yx Zx
  *          *  Xy Yy Zy
  * [x y z]     Xz Yz Zz]
  *
  * @param (const Vector3<TYPENAME>&) u - row vector
  * @param (const Matrix3<TYPENAME>&) m - matrix
  * @return Vector3<TYPENAME>
  */
 template <typename TYPENAME>
 inline Vector3<TYPENAME> operator*( const Vector3<TYPENAME>& u, const Matrix3<TYPENAME>& m )
 {
    return Vector3<TYPENAME>( (u.x*m.Xx) + (u.y*m.Xy) + (u.z*m.Xz), 
                              (u.x*m.Yx) + (u.y*m.Yy) + (u.z*m.Yz), 
                              (u.x*m.Zx) + (u.y*m.Zy) + (u.z*m.Zz) );
 }
 
 #endif // VECTOR3_H
 
 
 
 #endif // MATRIX_H